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Abstract 

A leading interpretation of the electron/positron excesses seen by PAMELA and ATIC 
is dark matter annihilation in the galactic halo. Depending on the annihilation channel, the 
electron/positron signal could be accompanied by a galactic gamma ray or neutrino flux, 
and the non-detection of such fluxes constrains the couplings and halo properties of dark 
matter. In this paper, we study the interplay of electron data with gamma ray and neutrino 
constraints in the context of cascade annihilation models, where dark matter annihilates 
into light degrees of freedom which in turn decay into leptons in one or more steps. Electron 
and muon cascades give a reasonable fit to the PAMELA and ATIC data. Compared to 
direct annihilation, cascade annihilations can soften gamma ray constraints from final state 
radiation by an order of magnitude. However, if dark matter annihilates primarily into 
muons, the neutrino constraints are robust regardless of the number of cascade decay steps. 
We also examine the electron data and gamma ray/neutrino constraints on the recently 
proposed "axion portal" scenario. 



1 Introduction 



Recent observations by PAMELA pQ and ATIC [2] strongly suggest a new primary source of 
galactic electrons and positrons. Three leading interpretations of the PAMELA/ATIC excesses 
are astrophysical sources [3J, decay of dark matter [I], and annihilation of dark matter [SI El d El 
|9]. While the current PAMELA/ATIC data cannot distinguish between these possibilities, one 
expects that the correct scenario will ultimately be determined with the help of complementary 
data from synchrotron, gamma ray, and neutrino telescopes, as well as collider and direct detection 
experiments. 

One piece of data that points toward an annihilation interpretation is the WMAP Haze |10j . 
an apparent excess of synchrotron radiation coming from the galactic center. Dark matter an- 
nihilation into charged particles is uniquely positioned to explain the Haze [TT1 H2] . If n is the 
dark matter number density near the galactic center, then the synchrotron signal for dark matter 
annihilation scales like n 2 , while the signal for dark matter decay scales only as n. (Astrophysical 
signals also roughly scale like n.) Given the normalization of the PAMELA/ATIC excess, the n 2 
scaling is favored to explain the size of the Haze anomaly [13] . 

On the other hand, the same n 2 versus n logic implies that the dark matter annihilation 
interpretation is more strongly constrained by the absence of gamma ray or neutrino excesses 
from the galactic center. While these constraints are dependent on the Milky Way dark matter 
halo profile, there are already strong bounds on the annihilation interpretation for strongly peaked 
halos [HJ ESI ESI EH EH]- Therefore, it is worth exploring dark matter annihilation scenarios in 
detail to understand how robust the tension is between explaining PAMELA/ATIC/Haze and 
satisfying other bounds. 

Given the absence of anti-proton [19] or gamma ray [201 [2U [22] excesses, the dark matter 
annihilation scenarios favored to explain PAMELA / ATIC involve annihilation into electrons and 
muons. However, dark matter need not annihilate into leptons directly. There are a variety of 
"cascade annihilation" models where dark matter annihilates into light resonances which in turn 
decay into electrons or muons. These light resonances can lead to nonperturbative enhancements 
[231 [21] of the dark matter annihilation rate in the galactic halo, providing the large boost factors 
necessary to explain PAMELA/ATIC [5j [6j [T| . Also, annihilation into light fields gives a kinematic 
explanation for why dark matter annihilation preferentially yields light leptons [121 El U\ ■ Previous 
studies of cascade annihilation models appear in Refs. [H ES]. 

In the present context, these cascade annihilation scenarios are interesting because they have 
the potential to explain PAMELA/ATIC while weakening constraints from gamma rays, as mea- 
sured by atmospheric Cerenkov telescopes like H.E.S.S. [2U1 Ell E2]- The reason is that gamma 
ray experiments are directly sensitive to the primary injection spectra, and cascade annihila- 
tions yield softer and smaller injection spectra of gamma rays from final state radiation (FSR). 
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PAMELA/ATIC sees electrons and positrons through the filter of charged cosmic ray trans- 
port, a process which introduces large uncertainties. Considering also the uncertainties in the 
highest energy ATIC data, we find that softer spectra of primary leptons can still explain the 
PAMELA/ATIC excesses. 

For cascade annihilations that terminate in muons, there is also an irreducible source of galactic 
neutrinos, which can be observed as an upward-going muon flux on earth, for example, by water 
Cerenkov detectors like Super-Kamiokande (Super-K) [25]. While cascades soften the neutrino 
spectrum, we will see that the final constraints from neutrinos are rather insensitive to the number 
of cascade steps, and may provide the most robust bound on muon cascade scenarios. 

The organization of this paper is as follows. In the next section, we define our framework 
for analyzing signals of dark matter through cascade annihilations, with details of the cascade 
energy spectra given in Appendix |A] In Section [3j we find the best fit dark matter masses and 
annihilation cross sections for various cascade scenarios given the PAMELA/ATIC data. We 
consider H.E.S.S. gamma ray bounds from FSR in Section H] and Super-K neutrino bounds in 
Section [5j In Section El we study a particular cascade annihilation scenario called the axion 
portal [7] , and present a less constrained "leptonic" version in Appendix Q31 Conclusions are 
given in Section [TJ 

2 Cascade Annihilations 

If dark matter is a thermal relic, then it will have at least one annihilation mode into standard 
model fields, since in the early universe, the dark matter annihilation channels keep dark matter 
in thermal equilibrium with the standard model until freezeout. However, dark matter need not 
annihilate into standard model particles directly; it can annihilate into new (unstable) resonances 
which in turn decay into standard model fields. As long as the new resonances are sufficiently 
broad, then dark matter will be in close enough thermal contact with the standard model for a 
freezeout calculation to be valid. 

These "cascade annihilations" can occur in one or more steps, as shown in Figure HJ and 
there are a variety of motivations for such cascade scenarios. Since the direct coupling between 
dark matter and the standard model can be small while still achieving the desired thermal relic 
abundance, such models can have reduced direct detection cross sections compared to generic 
weakly interacting massive particles (WIMPs) [2H]. Cascade annihilations can arise whenever 
new light resonances have large couplings to dark matter, such as in exciting dark matter |27j . 
In the context of PAMELA/ATIC, there has been interest in using light resonances to provide 
large enhancement to the galactic annihilation rate through the Sommerfeld effect or bound state 
(WIMPonium) formation [SI El El- These light resonances can also explain the lepton-richness of 
dark matter annihilation through kinematic thresholds [61 [7J . 
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Figure 1: For an n-step cascade annihilation, dark matter x annihilates into 4> n 4> n . The cascade 
annihilation then occurs through cf> i+ i — > (p^i (i — 1, ■ ■ ■ , n — 1), and in the last stage, (pi decays 
into standard model particles. The figure represents the cases where <pi — > £ + £~ . 

Cascade annihilations give softer primary spectra for the annihilation products than direct an- 
nihilation. As reviewed in Appendix |Aj for (scalar) cascade annihilations involving well-separated 
kinematic scales, the approximate primary spectra can be calculated using a simple convolu- 
tion formula. The energy spectra are conveniently parametrized in terms of the energy fraction 
x n = E n /m DM , where E n is the final state energy after the n-th step of the cascade, and m DM is 
the dark matter mass. If dN/dxo is the normalized primary spectrum from direct annihilation, 
then the normalized primary spectrum from a 1-step annihilation dN/dxi is 



This convolution formula can be iterated to obtain the normalized primary spectrum for an n-step 
cascade scenario. Eq. ([I]) is also reasonably representative of non-scalar cascades, since indirect 
detection measurements are unpolarized. 

While cascade annihilations give softer spectra, they typically yield a higher multiplicity of 
annihilation products. The final primary spectrum dN/ dx n must take into account the multiplic- 
ity of annihilation products per dark matter annihilation, and in a typical n-step cascade, the 
final state multiplicity scales like 2 n . 

An important exception, however, is in cases involving singularities, where cascade annihi- 
lations yield both a softer spectrum and a lower multiplicity. For example, FSR from charged 
leptons has a collinear singularity that is regulated by the lepton mass. Integrating over the 
singularity gives a photon spectrum dN 1 /dx that is proportional to lnQ/m^, where Q is the total 
energy available for radiation. For direct annihilation Q ~ 2mDM, while for cascade annihilations 
Q ~ nil, where mi is the mass of the resonance in the last stage of the annihilation <pi — > £ + £~~ . 
Light enough 0i fields can give a dramatic reduction in the In Q/rri£ factor and thus the FSR 
photon yield. 




(1) 
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In this paper, we only consider cascade annihilations that terminate in electrons or muons. 
The same analysis, however, could be repeated for charged pions or taus by changing the direct 
annihilation spectra. We expect the results for an n-step charged pion cascade to be similar to an 
(n + l)-step muon cascade. Cascades involving taus will face stronger gamma ray bounds because 
of the 0(1) fraction of 7r°s in tau decays. For simplicity, we only show plots for direct, 1-step, 
and 2-step cascades. 

3 PAMELA/ATIC Spectra 

The PAMELA satellite experiment |28j observed an anomalous source of galactic positrons in the 
energy range 10-100 GeV through a measurement of the positron fraction $ e +/($ e + + $ e -) PQ. 
The ATIC balloon experiment [29J is not capable of charge separation, but observed a peak in 
the total electron plus positron flux between 100-1000 GeV in a (<3> e + + $ e -) measurement [2]. 
Intriguingly, both excesses can be described by a single new source of galactic electrons and 
positrons, and here we study the goodness of fit for dark matter cascade annihilations. The 
primary electron/positron spectra for n-step electron and muon cascade annihilations are given 
in Appendices IA.1I and IA.2I Once the primary spectrum is known and a dark matter halo profile 
assumed, we can propagate the electrons and positrons through the Milky Way and compare with 
the PAMELA/ATIC data. 

We follow the analysis of Ref . [30] , which assumes that galactic electrons and positrons can be 
described by a diffusion-loss process. In the turbulent galactic magnetic fields, electrons/positrons 
diffuse within a fiducial region around the galactic disk and escape the galaxy outside that region. 
An energy loss term incorporates the physics of inverse Compton scattering (ICS) and synchrotron 
radiation. Dark matter annihilation is represented by a source term proportional to the square 
of the dark matter halo density. Since the energy loss time is much shorter than the age of the 
galaxy, the electron/positron system is assumed to be in steady state. 

Taking ip e - (x, E) to be the galactic electron number density per unit energy, the diffusion-loss 
equation is 

K e 5 V 2 W (x , E) + -^(^- (x, E)^j + q(x, E) = 0, (2) 

where e = E/GeV, K and S parametrize the (energy dependent) diffusion, te is a characteristic 
energy loss time, and q(x, E) is the electron source term for dark matter annihilations. The 
same equation also holds for the positron number density per unit energy ip e +(x,E). The elec- 
tron/positron densities ip e -/ e + (x, E) are assumed to have vanishing boundary conditions on the 
surface of a cylinder of height 2L and radius R. We consider the three benchmark models from 
Ref. [30J , which are summarized in Table [TJ 
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R (kpc) 


L (kpc) 


K (kpc 2 /Myr) 


5 


t e (sec) 


MED 


20 


4 


0.0112 


0.70 


10 16 


Ml 


20 


15 


0.0765 


0.46 


10 16 


M2 


20 


1 


0.00595 


0.55 


10 16 



Table 1: Diffusion-loss parameters for the three benchmark models (MED, Ml, and M2) for 
electron/positron propagation. 

The electron/positron source term is given by 

where p(x) is an assumed dark matter halo profile, moM is the dark matter mass, (av ) is the 
average dark matter annihilation cross section in the galactic halo, and dN e /dE is the electron 
energy spectrum per dark matter annihilation, rj — 1 if dark matter is self-conjugate (e.g. a 
Majorana fermion), while 77 = 2 if not (e.g. a Dirac fermion). We consider three spherically 
symmetric benchmark halo profiles (cored isothermal [21], NFW [32], and Einasto [33]) with 
r = 8.5 kpc and p = 0.3 GeVcm -3 : 

\ _|_ (r If )2 

P(r)lsothermal = P® ~. , 7 ~, , T c = 5 kpc, (4) 

1 + {r/r c ) 2 

p(r) N FW = Pq— ( 1 1 +r °/ r ) , r c = 20 kpc, (5) 
r \ 1 + r/r c J 

P(0Einasto = PQ eX P ( ~ ) - f — 

^-body simulations suggest that Einasto and NFW are more realistic profiles for r 1 kpc. 
Within the inner region, however, there is considerable uncertainty, and we include the cored 
isothermal profile to explore the possibility of a less peaked distribution. 

Once the source term is specified, Eq. (T2]) can be solved using the methods of [301 El], and the 
electron/positron intensities (fluxes per energy per solid angle) at the earth due to dark matter 
annihilations are given by 

• (DM) 

-(E) = ^p^^ e - /e+ (x & ,E), (7) 



a = 0.17, r c = 20 kpc. (6) 



Uj ^e~/e+ , -D e) astro 



where x & is the location of the solar system and -B e ,astro is an astrophysical boost factor. Since 
the annihilation rate is proportional to the squared density of dark matter P5 M , clumping of 
dark matter tends to increase the local electron/positron flux, and the -B e ,astro factor accounts for 
differences from the assumed smooth halo profile, as well as uncertainties in p Q . Of course, in the 
presence of dumpiness, the true electron/positron spectrum is also modified [35] . 
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Galactic cosmic rays are a known source of background primary electrons. Background sec- 
ondary electrons and positrons arise, e.g., from collisions of cosmic ray protons with interstellar 
gas. In principle, the spectra of background electrons/positrons are correlated with the diffusion- 
loss parameters for electrons/positrons, but for simplicity we will use the parameterization of 
background primaries and secondaries from Ref. [36J: 

GeV cm s sr , (8) 



dEdn l + lle°- 9 + 3.2e 2 - 15 

d^ c) 0.70 e - 7 



GeV^cm-'s-'sr- 1 , (9) 



dEdtt 1 + 110£ L5 + 600£ 2 - 9 + 580£ 4 - 2 

GeV^cm^s^sr- 1 , (10) 



dEdVL 1 + 650 e 2 - 3 + 1500 e A - 2 

where again e = E/GeV. To treat background uncertainties, we will marginalize over the nor- 
malization and overall slope of the background in our analysis: 

$ (back. fit) = A _ £ P_ ^(prim) + $ (sec)^ ^ $ (back. fit) = ($(^ ^ ( n ) 

where we allow < A± < oo, -0.05 < P± < 0.05, as in [5]. 

In the limit that Eq. ([I]) holds, any given dark matter cascade topology has just two free 
parameters: the dark matter mass m^M and the annihilation cross section (av). Following the 
literature, we normalize the cross section to the value that leads to the right relic thermal abun- 
dance (av)o ~ 3 x 10~ 26 r] cm 3 s" 1 and express our results in terms of an effective boost factor 

B = -B e , astro 7 T~, (12) 

which includes both the deviation from the naive thermal freezeout cross section and dark matter 
clumping. Using only the statistical error bars, we perform a chi-squared fit of the derived 
electron/positron intensities to the PAMELA e + /(e + + e~) and ATIC e + + e~ data, treating todm 
and B as free parameters. We use only E ^ 10 GeV bins for the PAMELA data, as the lower 
energy bins are strongly affected by solar modulation effects (see the discussion in [1]). 

The results of the fit for direct, 1-step, and 2-step annihilations are shown in the case of 
electron final states in Figure [2] and muon final states in Figure |3j (The results do not depend 
on whether the dark matter particle is self-conjugate or not.) Each plot corresponds to a definite 
cascade annihilation pattern with a definite halo profile, and shows la and 2a contours for three 
propagation models in the mjyu-B plane. The best fit values for mDM and B, as well as the \ 2 
values, are also given, where x 2 follows the chi-squared distribution with 22 degrees of freedom 
(7 PAMELA + 21 ATIC — 6 fitting parameters: ttidm, B, A±, P±). To evaluate goodness of fit, we 
also give p- values in Table [2, where we have chosen the propagation model giving smallest x 2 f° r 
each plot of Figures [2] and O As we can see, the fits are reasonable for all the cases presented. 
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Figure 2: The best fit regions for the dark matter mass ttl-qm and boost factor B in the cases 
of direct, 1-step, and 2-step annihilations into e + e~ for different halo profiles and propagation 
models. The best fit values are indicated by the crosses, and the contours are for la and 2a. 
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Figure 3: The same as Figure [2] but for annihilations into . 
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Table 2: The p- value for the best propagation model for each plot in Figures [2] and El 
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Figure 4: The predicted e ± intensities compared to the PAMELA (left) and ATIC (right) data 
for direct (solid), 1-step (dashed), and 2-step (dotted) annihilations into electron final states. 
The NFW halo profile and the MED propagation model are chosen, and the e ± backgrounds are 
marginalized as described in Eq. ffTT]) . Note that we fit the PAMELA data only for E ^ 10 GeV 
because solar modulation effects are important at lower energies. 
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Figure 5: The same as Figure H] but for annihilations into muon final states. 
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In Figures H] and we show the comparisons of the e ± intensities with the PAMELA and 
ATIC data in the cases of electron and muon final states. Solid, dashed, and dotted lines represent 
direct, 1-step, and 2-step annihilations, and we have chosen the NFW halo profile and the MED 
propagation model for illustrative purposes. We clearly see the trend that more steps in cascades 
lead to flatter e ± spectra, but that the fits are good for all cases shown. Note that we have not 
optimized the propagation model here. Adjusting the propagation model can lead to a better fit 
in certain cases, especially the 2-step annihilation into 

Summarizing the analysis of the e ± data, we find: 

• The PAMELA and ATIC data require the dark matter mass and boost factor in the region 
jtidm = O(TeV) and B = 0(1000), which is consistent with earlier analyses on the direct [5] 
and 1-step [8] cases. More steps in the cascade lead to larger values of mDM and B, and 
roughly speaking, both B and m^M scale as 2 n . The reason is that the peak location in 
the ATIC data sets the scale for the final state e energy, and since the average e energy 
over m DM in an n-step cascade scales like 1/2™, 7«dm must increase by 2™ to keep the peak 
location fixed. Similarly, the annihilation signal scales like BN e /m\ iM (assuming fixed Pdm); 
where N e = 2 n+1 is the final state e ± multiplicity. Thus, to keep the PAMELA/ATIC rate 
fixed, B must scale like 2™. 

• The fits do not become much worse by going to multiple steps, due to uncertainties in the 
highest energy ATIC data and uncertainties in the e propagation model. In particular, 
2-step annihilations still fit the data reasonably well for both electron and muon final states. 
The required boost factors are rather large in the case of muon final states: B of a few 
thousand. Such large boost factors may come from both astrophysics, e.g. uncertainties in 
Pq and nearby clumps of dark matter, and particle physics. 

• Uncertainties in halo profiles and propagation models do not significantly affect the dark 
matter mass and the boost factor. Errors from these uncertainties are mostly of 0(10%) 
and at most a factor of 2. The 2a ranges of the fits then determine m^M and B up to a 
factor of a few. 

Recent measurements at H.E.S.S. [37] of the electron plus positron flux above 600 GeV [38] are 
qualitatively consistent with the ATIC spectrum. Given the large systematic energy uncertainties 
and hadronic background, we do not use the H.E.S.S. data in our fits, although we remark that 
the observed steepening of the spectrum places some bounds on very long cascade decays and 
may disfavor spectra with a hard cutoff such as direct e annihilation. 
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4 Gamma Ray Constraints 



When dark matter annihilates into charged leptons, there is a primary source of gamma rays 
coming from FSR. Various gamma ray telescopes have looked for excess gamma rays coming from 
the galactic center, and the null result of such searches puts bounds on dark matter annihilation 
into charged leptons. An additional effect that is beyond the scope of this paper is ICS, where 
electrons/positrons from dark matter annihilation lose energy by upscattering starlight photons 
into gamma rays. The rate of ICS photon production depends on the modeling of galactic 
starlight, and we here focus only on the bounds from FSR. For an early analysis of FSR in dark 
matter annihilation, see [39] . 

There is negligible energy loss as gamma rays propagate from the galactic center to the earth. 
The total power (flux per energy) on earth depends on the dark matter halo profile through 

= Castro {av) jAn ^ (13) 

dE 87r7/mjV 1 dE ' v ' 

where I? 7i astro is an astrophysical boost factor for photons that may differ from -Bistro, and 
dN^/dE is the photon energy spectrum per dark matter annihilation. The energy spectrum 
from FSR is reviewed in Appendix IA.31 and we also include the effect of radiative muon decays 
as described in Appendix IA.41 AQ is the solid angle integration region, and J is the average 
line-of-sight-integrated squared dark matter density for a given halo model 

J = -^f dnf ds P {xf. (14) 

J An Jline-of-sight 

The strongest bounds on FSR gamma rays come from atmospheric Cerenkov telescopes, but 
the way these experiments extract gamma ray signals affects the final dark matter annihilation 
bounds. To enable background subtraction, these telescopes operate either in on-off mode or 
wobble mode, meaning the effective J exposure is [40] 



Jeff Jon— source Joff— source- (1^) 

By definition, J er r is smaller than J on -source, and neglecting the Jog-source contribution gives bounds 
that are too aggressive. This means that one must know both the on- and off-source integration 
regions to derive a bound on the cross section. For shallow dark matter halo profiles, there can 
be large cancellations in the value of J en -, and in principle a stronger bound could be obtained 
using the raw unsubtracted data. 

We set bounds on FSR using three H.E.S.S. gamma ray data sets. The first two are obser- 
vations of the Galactic Center (GC) [20] and the Galactic Ridge (GR) [21]. Neither is ideal for 
dark matter observations because of the large contamination from gamma ray point sources and 
molecular gas, and in principle one should put bounds on a dark matter signal after subtracting 
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GC (on) 


GC (off) 


GC (eff) 


GR (on) 


GR (off) 


GR (eff) 


Isothermal 


10 


10 


0.028 


10 


10 


0.019 


NFW 


1.1 ■ 10 4 


3.6- 10 2 


1.1 ■ 10 4 


1.8 • 10 3 


4.3 ■ 10 2 


1.4 • 10 3 


Einasto 


5.8- 10 3 


7.3 ■ 10 2 


5.1 • 10 3 


2.3 • 10 3 


8.7- 10 2 


1.5 ■ 10 3 



Table 3: J values for GC and GR gamma ray observations (on-source, off-source, and effective) 
in units of GeV 2 cm -6 kpc. 

both these foregrounds. Since such subtractions are not available, we derive conservative bounds 
by insisting that the dark matter signal does not exceed any of the H.E.S.S. data points by more 
than 2a. For the GC and GR samples, values of J for the three dark matter halos in Eqs. (jH El ED 
are shown in Table [3] for both the on-source and off-source regions. 

The GC data set comes from the inner 0.1° of the galaxy with a solid angle integration of 
AQ = 1 x 10~ 5 , corresponding to the gamma ray source HESS J1745-290. This sample was taken 
in wobble mode, and the off-source region corresponds to a ring at a distance of 1.4° from the 
GC. Apart from the off-source subtraction, no other corrections were made to the data, so the 
data points include both the HESS J1745-290 point source as well as any putative dark matter 
signal. 

The GR data set comes from the region \£\ < 0.8°, |6| < 0.3° in galactic coordinates, with 
foreground point sources HESS J1745-290 and GO. 9+0.1 subtracted. The GR sample was taken in 
on-off mode, and the region \£\ < 0.8°, 0.8° < |fe| < 1.5° was used for background subtraction. In 
the GR data, H.E.S.S. finds that the gamma ray emissions are spatially correlated with molecular 
gas traced by CS emission lines, but a molecular gas foreground component is not subtracted 
from the data. The GR bounds are expected to strengthen if one were to subtract a molecular 
gas component. 

H.E.S.S. looked more directly for dark matter annihilation in the Sagittarius dwarf spheroidal 
galaxy (Sgr dSph) [22]. Sgr dSph has negligible foregrounds and is thought to be dark matter 
dominated. A 95% C.L. model-independent bound on the gamma ray flux was obtained for 
£ 7 > 250 GeV: 

$ 7 (£ 7 > 250 GeV) < 3.6 x 10~ 12 cm" 2 s" 1 , (16) 

with a solid angle integration of AQ = 2 x 10~ 5 . Since the Sgr dSph data was taken in wobble mode 
and the contribution from the Sgr dSph halo is negligible in a 1.4° ring, we use J e s = ^on-source- 
The value of J e e strongly depends on the halo profile of Sgr dSph. For example, an NFW profile, 
a large core profile, and a small core profile quoted in [15] lead to 

rSgrdSph _ 7 o iq2 rSgrdSph _ , , ]() 2 rSgrdSph _ p a -i r,4 /-irA 

J NFW — /.O X 1U , ^LargeCore " 1J x iu I J SmallCore — X iU > \ lt ) 

in units of GeV 2 cm -6 kpc. To derive bounds, we consider NFW and large core profiles. 
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Figure 6: Constraints from gamma ray observations, GC (solid), GR (dashed), and Sgr dSph 
(dotted), in the ui-dm-B plane for direct, 1-step, and 2-step annihilations into electron final states. 
All the constraints, as well as the best fit region for PAMELA/ATIC (MED propagation), are 
plotted assuming -B e ,astro = -B 7 astro . For cascade annihilations, each of the GC, GR, and Sgr dSph 
constraints consist of two curves, with the upper (blue) and lower (red) curves corresponding to 
mi = 100 MeV and 1 GeV, respectively. Note that the constraint lines in the cored isothermal 
case are above the plot region, and that the halo profiles for Sgr dSph are given in the legends. 
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Figure 7: The same as Figure [6] but for muon final states. Also included are constraints from 
neutrino observations (dot-dashed) assuming B VjastIO = -B e ,astro- For cascade annihilations, the 
upper (blue) and lower (red) curves now correspond to mi = 600 MeV and 1 GeV, respectively. 
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In Figures and [7J we show the resulting constraints from the GC, GR, and Sgr dSph gamma 
ray observations in the 7Bdm--B plane. We also superimpose the la and 2a contours reproducing 
the PAMELA/ATIC data for the MED propagation model; see Figures [2] and [31 In order to plot 
the Sgr dSph and Milky Way bounds on the same plane, we associate the Sgr dSph small core 
profile with the Milky Way cored isothermal profile, and the Sgr dSph NFW profile with the 
Milky Way NFW and Einasto profiles. The ragged lines in the GC and GR constraints come 
from the binning of the H.E.S.S. data. 

We note that here we have drawn the gamma ray constraints and e contours assuming a 
common -B 7jas tro and -Bistro- This is most likely not the case. For example, if dark matter 
clumping decreases toward the galactic center, then the gamma ray flux will decrease compared 
to the local positron flux. If the astrophysical boost factor for e is larger than that for 7, then 
the gamma ray constraints become weaker by a factor of B etastro / '-Bistro compared to those shown 
in Figures [6] and [7J 

The bounds from Sgr dSph can be modified if there are nonperturbative enhancements to the 
dark matter annihilation cross section. Since the velocity dispersion of dark matter in Sgr dSph is 
"^Sgr ~ 10 km/s [UJ, as opposed to v MW ~ 200 km/s in the Milky Way, the boost factor relevant 
for Sgr dSph may be larger than that for electrons/positrons if part of the boost factor arises 
from the Sommerfeld or bound state enhancement. This would make the Sgr dSph constraint 
stronger than what is naively read from Figures [6] and UJ 

Summarizing the analysis of the gamma ray constraints, we find: 

• For FSR gamma rays from the galactic center region, the GR data gives somewhat stronger 
constraints than the GC data. For direct annihilation, this disfavors the NFW and Einasto 
profiles, as also seen in Refs. [HJ [T3] . The Sgr dSph data is not constraining unless we were 
to take a highly peaked halo such as the small core profile. 

• The constraints from FSR photons are weaker in the cascade annihilation case than in the 
direct annihilation case. This is because cascade annihilations give smaller photon yield 
at high energies, as discussed in Appendix IA.31 The smaller the mass is, the weaker 
the constraints become. The constraints, however, do not become weaker by increasing the 
number of steps, as can be seen by comparing the best fit region for the e data with the 
gamma ray constraints. This is because the reduction in the photon yield depends only on 
the 0i mass and not on the number of cascade steps, and the softening of the gamma ray 
spectra is compensated by the increase in the best fit dark matter mass and boost factor. 

• The constraints are very weak for shallow halo profiles such as the cored isothermal profile 
for the Milky Way in Eq. (jl]). This is particularly true for the GC and GR data because 
the background subtraction due to on-off or wobble mode operation also subtracts (most of) 
the signals from dark matter annihilations. The constraints from these data, therefore, are 
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rather weak as long as the halo profile is relatively flat within about 100 pc of the galactic 
center. A better bound may be obtained if we could use the unsubtracted data. 

While we have focused only on FSR photons in our analysis, we wish to make a few comments 
about ICS, the WMAP Haze, and radio bounds. In the context of dark matter, the WMAP 
Haze arises because electrons from dark matter annihilation emit synchrotron radiation in the 
galactic magnetic fields. The total synchrotron power — and hence the size of the predicted WMAP 
Haze signal — depends on whether these electrons can lose energy via non-synchrotron channels. 
Ref. [8] found that in order to be consistent with the large boost factors necessary to explain 
PAMELA/ATIC, one had to assume a larger rate for ICS compared to earlier Haze analyses |llj . 
Given the uncertainty in galactic starlight and the dark matter halo profile, it is consistent 
to conservatively ignore potential bounds from ICS, but since the WMAP Haze is one of the 
motivations for considering dark matter annihilation, strictly speaking one should verify that the 
assumed electron energy loss mechanisms can yield the WMAP Haze while satisfying ICS photon 
bounds. 

That said, we do not expect much variation in the ICS bounds between direct and cascade 
annihilation scenarios. ICS is calculated from a steady state configuration of charged particles, so 
to the extent that the PAMELA/ATIC data already normalizes the steady state electron/positron 
densities, the ICS yield should be similar regardless of the annihilation scenario. Therefore, the 
recent analysis of Ref. [8] should be representative of generic multi-step cascade annihilation 
scenarios. This is similar in spirit to the WMAP Haze, in that the Haze requires a source of 
charged particles to generate the synchrotron signal, but the precise particle energy distribution 
has only a secondary effect. 

Finally, there has also been recent suggestions [T5J [T6] of a possible tension between a dark 
matter annihilation interpretation of the WMAP Haze and bounds from 408 MHz radio obser- 
vations in the inner 4" of the Milky Way [42]. Such bounds assume that the steep halo profiles 
necessary to generate the Haze at a latitudinal distance between 5° and 30° can be extrapolated 
to sub-parsec distances away from the galactic center. There are a number of reasons to distrust 
such an extreme extrapolation of the dark matter halo, including possible effects of baryons [43] 
and hierarchical mergers [44] to soften cuspy behavior. At minimum, iV-body simulations [43] do 
not have resolution to such small scales. Therefore, we find no reason to disfavor a dark matter 
annihilation scenario on the basis of the 408 MHz radio bound. Note that the analysis of Ref. [45] 
using an all-sky radio model [47] finds only relatively mild synchrotron constraints on TeV-scale 
dark matter after masking the inner 15° x 15° of the galaxy. 
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5 Neutrino Constraints 



When dark matter annihilates into muons, there is an irreducible source of neutrinos. Neutrinos 
produced in the galactic center oscillate as they travel towards earth, and if they are muon-type 
neutrinos when they collide with rock in the earth's crust, they can create an upward-going flux 
of muons. These muons could be observed by water Cerenkov detectors, and the absence of such 
observations puts bounds on the dark matter annihilation rate into muons. For muon cascades, 
there is no high energy neutrino source from dark matter that accretes in the sun and earth, 
because the muons from dark matter annihilation are stopped before they decay [38] . 

Since neutrinos have negligible energy losses as they traverse the galaxy, the muon-neutrino 
flux incident on earth is 

%^ = E^^f-M^An^, (is) 

where i runs over the neutrino flavors, B vastvo is the astrophysical boost factor for the neutrino 
signal which could differ from -Bistro, and P ViV is the probability that V{ has oscillated into 



P v „ = 0.39, i> = 0.21. (19) 



There is an analogous formula for z/ M , but in dark matter annihilations the and fluxes are 
equal. The primary neutrino spectra are given in Appendix IA.2I 

We now calculate the resulting upward-going muon flux following the analysis of Ref. [50] . 
An incident neutrino of energy E v creates muons of energy E^ according to the neutrino-nucleon 
scattering cross sections (Jun^^x- For the propagation of created muons, we use an approximate 
energy loss parameterization 

-£ = p m *(-a-/3E), (20) 

with "standard rock" values a = 2 x 10 -6 TeVcm 2 /g and /3 = 4 x 10~ 6 cm 2 /g; p mat will cancel in 
the final muon flux expression. In this approximation, a muon of starting energy can travel a 
distance 

L(E„ E thTCS .) = - J- In f °[ ) (21) 

before its energy drops below the muon detection threshold E thTes . 
For a given d$„ JdE v , the observed muon flux is 

rm-Du fm D M I d<& v do v N-,ii~X d<&v d(J 9 N-^n+x\ , , „ 

dEj dE v n N ^ + 'fj^T L(E„ E thies .), (22) 

JEthres. Je» \ dE v dE^ dE v dE^ J 

where = Pmat/^iv is the nucleon number density in the earth's crust, and N refers to an 
average nucleon. We calculate the neutrino-nucleon scattering cross sections assuming an equal 
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3° 5° 10° 15° 20° 25° 30' 



Isothermal 
NFW 
Einasto 



10 10 9.7 9.0 8.2 7.4 6.6 
340 190 84 51 35 26 20 
640 370 160 91 60 43 32 



Super-K (95% C.L.) 



2.70 4.82 6.43 10.6 11.2 17.6 19.5 



Table 4: J values for neutrino observations in units of GeV 2 cm 6 kpc, and Super-K 95% C.L. 
flux limits in units of 10~ 15 cm~ 2 s _1 . 

fraction of protons and neutrons in rock, using CTEQ5M parton distribution functions [51] to 
include the effect of sea quarks, and retaining the full W boson propagator in the cross section. 

Super-K [52J placed 90% confidence bounds on the upward-going muon flux [25] in various 
cone sizes ranging from 3° to 30° around the galactic center, with -E t hres. = 1.6 GeV. To be more 
conservative, we consider 95% confidence bounds |53j, as shown in Table HI The relevant values 
of J for the three dark matter halos in Eqs. (JU, El ED are a ^ so shown in Table 0} To derive a bound 
on the annihilation rate, we insist that the predicted flux does not exceed the 95% confidence 
bound for any of the Super-K cone sizes. The neutrino constraints for muon cascade scenarios 
appear in Figure 0, assuming a common boost factor for the electron, gamma ray, and neutrino 
signals. 

As observed in Ref. [48], for sufficiently small dark matter masses, the observed muon flux 
is nearly independent of the dark matter mass. The reason is that both the neutrino-nucleon 
scattering cross section and the muon range scale like energy, but the dark matter number density 
squared (and hence the annihilation signal) scales like I/ttidmj so the final observed flux is simply 
related to the normalized second-moment of the neutrino energy spectrum. Hence the exclusion 
limit for the neutrino boost factor is approximately flat in todm for small enough dark matter 
masses. As the dark matter mass increases, the average neutrino gets harder, and the neutrino- 
nucleon cross section grows less steeply because of the W boson propagator. In addition, the 
energy-dependent term in dE/dL begins to take effect, relaxing the neutrino bounds for high 
dark matter masses. 

Summarizing the analysis of the neutrino constraints, we find: 

• Since the neutrinos are softer in cascade annihilations, the bounds on the boost factor 
are weaker than for direct annihilation. However, the electrons are also softer in cascade 
annihilations, so the PAMELA/ATIC best fit mass and boost factor rise. Put together, the 
neutrino tension increases marginally as the number of cascade steps increase. For direct 
annihilation, our results agree qualitatively with [T71 \TE\ . 

• Super-K considered solid angles as large as 30°, so a large fraction (~ 30% to ~ 70%) of 
the total dark matter annihilation signal is contained within the observed region. While 
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Figure 8: In the axion portal, fermionic dark matter annihilates dominantly into a scalar s and 
a pseudoscalar "axion" a. The scalar then decays as s — > aa, and the axion decays as a — > £ + £~ . 
In the minimal axion portal, the axion dominantly decays into muons, but in the leptonic axion 
portal it can dominantly decay into electrons. These models are partway between a 1-step and a 
2-step cascade annihilation scenario. 

the Einasto profile is less peaked than NFW toward the galactic center, the Einasto bound 
happens to be stronger because of the large integration region. For the NFW and Einasto 
profiles, the dominant bounds come from the 10° cone, while for the cored isothermal profile 
they come from the 30° cone. 

• Since the gross structure of the dark matter halo is better understood than the halo density 
at the galactic center, the Super-K neutrino constraints are in some sense more robust 
than the H.E.S.S. gamma ray constraints. As in the case of gamma rays, the bounds are 
rather weak for the cored isothermal profile. They are, however, significantly stronger in the 
more realistic NFW and Einasto profiles. For cascade decays, the neutrino constraints are 
comparable to or stronger than the FSR gamma ray constraints, and highlight the tension 
in muon cascade scenarios. As with the gamma ray bounds, however, differences between 
the astrophysical boost factors -Bistro and B u ^ stro could alleviate the tension. 

6 The Axion Portal 

One well- motivated example of a cascade annihilation scenario is the axion portal [7J. In this 
scenario, dark matter is a TeV-scale particle that obtains a mass from spontaneous symmetry 
breaking. The spontaneous breaking of U(l)x yields a pseudoscalar "axion" a and a scalar 
"Higgs" s, and for fermionic dark matter the dominant annihilation channel is 

XX -> sa (23) 

(xx — ^ sa if x is a Majorana fermion). The scalar s dominantly decays as s — > aa, and if standard 
model leptons carry axial U(l)x charges, then a — ► £ + £~. Since a is a pseudoscalar, helicity 
suppression implies a will decay into the heaviest kinematically allowed lepton, which we assume 
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Minimal Axion Portal 



Leptonic Axion Portal 



Isothermal 
NFW 
Einasto 



7.9% 
7.9% 
8.3% 



9.2% 
9.2% 
8.8% 



Table 5: The p-values for the best propagation model for each plot in Figure [9] 

is either an electron or muon. An exchange of s can also provide the necessary enhancement 
of the annihilation cross section through the Sommerfeld and/or bound state effect. As shown 
in Figure [8], the axion portal effectively gives a one-and-a-half step cascade annihilation in the 
language used here. 

The simplest model for the axion portal — the minimal axion portal — is obtained if we iden- 
tify U(l)x with a Peccei-Quinn symmetry rotating two Higgs doublets. In this case, a has 
large hadronic couplings, and there are strong constraints on the axion mass from beam dump 
experiments and rare meson decays. Ref. [7J found that the preferred axion mass range was 
360 MeV ^ m a ^ 800 MeV, in which case a preferentially decays into muons. In the analysis 
here, we also consider a variant of the axion portal — the leptonic axion portal — where only lep- 
tons are charged under U(l)x- This possibility was mentioned in [7J, and is described in more 
detail in Appendix [Bj Dark matter in this model annihilates through the axion at associated 
with the leptonic symmetry, which does not have a coupling to hadrons. The absence of hadronic 
couplings allows the parameter range 2m e < m ae < 2m M , so at can preferentially decay into 
electrons. Note that the nonperturbative enhancement of the halo cross section is caused by the 
exchange of another scalar st and not by the axion at, so the bound of Ref. [51] does not exclude 
at masses smaller than ^100 MeV. 

In Figure [91 we show the best fit values for 771dm and B for the PAMELA/ATIC data for the 
three different diffusion models and three different halo profiles. In Figure [TJjl we compare the 
best fit regions to FSR and neutrino constraints. Here we assume that B is common for electron, 
gamma ray, and neutrino signals, so that the same qualifications as Figures El and [7J apply. To 
account for the fact that smaller masses for at are allowed, we are considering smaller values of 
m ai in Figure [10] than those of mi in Figure [6] For completeness, we also show the best fit spectra 
to the PAMELA and ATIC data in Figure [11], and the p- values of the fit in Table [5] As expected, 
the best fit values and qualitative features of the plots are partway between a 1-step and a 2-step 
cascade scenario. 

Depending on the axion mass, there are also potential gamma ray constraints from rare a — > 77 
and a — > 7r + 7r~7r° decay modes, and we can place bounds on the branching fractions to these 
modes. The a — > 77 mode is potentially dangerous for the entire axion mass range, but the 
a — > 7r + 7r~7r° mode only when m a > 3m n . The gamma ray spectra for these rare decays are given 



20 



IO 4 



B 



IO 2 



Minimal Axion Portal (Isothermal) 



Minimal Axion Portal (NFW) 



MED 3310 3850 34.0 
Ml 3310 3550 35.8 
M2 2470 3460 31.9 



10" 








B 


X 2 


MED 


3320 


3690 


34.7 


Ml 


4050 


4010 


46.7 


M2 


2460 


3560 


31.9 



10 3 




10 2 



10 3 



io 2 

10 4 to 2 



10 4 



10' 



Minimal Axion Portal (Einasto) 



MED 3310 3670 35.5 
Ml 4560 4220 58.7 
M2 2460 3480 31.7 




10 3 



10 2 h 
10 4 10 2 



10 ! 



10 4 



Leptonic Axion Portal (Isothermal) 



Leptonic Axion Portal (NFW) 



Leptonic Axion Portal (Einasto) 





m DM 


li 


x 1 


MED 


1340 


522 


31.2 


Ml 


1330 


478 


31.6 


M2 


1000 


429 


36.3 



B io 3 



to 2 



10 2 



10' 
m DM 



10 4 





m DM 


B 


x 1 


MED 


1350 


503 


31.2 


Ml 


1520 


478 


39.2 


M2 


1010 


458 


35.8 



10 3 




io 2 h, 

10 4 10 2 



1 o 4 



10' 



"DM 



B x 2 

MED 1320 481 31.4 

Ml 1750 538 50.7 

M2 1010 445 35.9 




10 3 
'"DM 



10 2 h 
10 4 10 2 



10 3 



10 J 



Figure 9: The best fit regions for the dark matter mass tbdm and boost factor B in the minimal 
axion portal (top row, a — > and leptonic axion portal (bottom row, — > e + e~) for different 

halo profiles and propagation models. The best fit values are indicated by the crosses, and the 
contours are for la and 2a. 



21 



Minimal Axion Portal (Isothermal) Minimal Axion Portal (NFW) Minimal Axion Portal (Einasto) 




10- 10 3 10" 10 2 10 3 10 4 10 2 10 3 10 4 



Leptonic Axion Portal (Isothermal) Leptonic Axion Portal (NFW) Leptonic Axion Portal (Einasto) 




m DM m DM m DM 



Figure 10: Constraints from gamma ray, GC (solid), GR (dashed), and Sgr dSph (dotted), and 
neutrino (dot-dashed) observations in the 7tidm--B plane in the minimal axion portal (top row, 
a — ► /U + /i~) and leptonic axion portal (bottom row, ai — > e + e~). All the constraints, as well 
as the best fit region for PAMELA/ATIC (MED propagation), are plotted assuming that B is 
common. Each of the GC, GR, and Sgr dSph constraints consist of two curves. For the minimal 
axion portal, the upper (blue) curve is m a = 600 MeV and the lower (red) curve is m a = 1 GeV. 
For the leptonic axion portal, the upper (blue) curve is m ai = 10 MeV and the lower (red) curve 
is m ae = 100 MeV, which differs from the choice in Figure [61 Note that the constraint lines in 
the cored isothermal case are above the plot region, and that the halo profiles for Sgr dSph are 
given in the legends. 
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Figure 11: The predicted e ± intensities compared to the PAMELA (left) and ATIC (right) data 
for the minimal axion portal (solid, a — > yU + /i~) and leptonic axion portal (dashed, ag — > e + e~). 
The NFW halo profile and the MED propagation model are chosen, and the e ± backgrounds are 
marginalized as described in Eq. (jlip . Note that we fit the PAMELA data only for E ^ 10 GeV 
because solar modulation effects are important at lower energies. 

in Appendix IA.5I For each halo profile and for the three H.E.S.S. data sets, we turn off FSR 
and find the values of the branching ratios where the H.E.S.S. bound is saturated for the best fit 
values of moM and B and the best propagation model (assuming B is common for e and gamma 
rays) . 

The branching ratio bounds are shown in Tables [6] and [7J The constraints from Sgr dSph 
match the estimates in [7j, but the GC and GR data imply a factor of 5-10 stronger bound. In 
the minimal axion portal, Br(a — > 77) ~ 10~ 5 -10~ 4 , which is safe from the bounds. There is 
somewhat more tension for a — > 7r + 7r~7r° where the estimated branching fraction was ~ 10~ 2 . 
Note, however, that the constraints from the GC and GR data are weak for any halo profile 
relatively flat within about 100 pc of the galactic center. In particular, there is no meaningful 
branching ratio constraint for the cored isothermal profile in Eq. (j4j). 

Gamma ray and radio bounds on the axion portal were also considered in Ref. [16], where 
potential s — » bb and s — > t + t~ decay modes were included. Since these additional s decay modes 
depend on details of the model, we do not consider them here. For FSR gamma rays, our results 
and conclusions agree qualitatively with [TB] . As mentioned already, we do not consider 408 MHz 
radio observations to place any meaningful constraint on the axion portal. 

7 Conclusions 

The possibility of indirect detection of dark matter has been considered for over 25 years |55j . 
but the annihilation rates expected from WIMP thermal relics are typically too small to give 
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Brfo, — ► / vo /N ) 




GC 


GR 


Sst dSnh 


NFW Best 
Einasto Best 


2.2 • 1(T 3 
4.8 • 1(T 3 


1.1 ■ 10~ 3 

1.2 ■ 10" 3 


2.5 ■ 10^2 


Br(a — > 7T + 7T~7T° 


) 


GC 


GR 


Sgr dSph 


NFW Best 
Einasto Best 


1.3-10" 
3.0- 10" 


" 2 7.5-10- 
- 2 7.9 ■ 10- 


3 5.2 • 10~ 2 

-3 



Table 6: Bounds from gamma rays on the branching fractions of a — > 77 and a — > 7r + 7r-7T° in the 
minimal axion portal (a — > These are obtained neglecting all other sources of gamma rays 

and correspond to the best fit values for m DM and £> and the propagation model giving smallest 
X 2 - The bounds assume an equal boost factor for e ± and gamma rays, and should be multiplied 
by -B e ,astro/-B7,astro if the boost factors differ. 



Br(a e -> 77) 


GC 


GR 


Sgr dSph 


NFW Best 


8.2 • lO" 3 


4.6 • 10" 3 


6.5 • 10-2 


Einasto Best 


1.8 • 10-2 


4.9 • 10- 3 





Table 7: The same as Table El but for the leptonic axion portal (an — > e + e ). The a — > n + ii n° 
mode is irrelevant in this case. 

appreciable gamma ray or neutrino fluxes from the galactic center unless a very peaked dark 
matter halo profile is assumed. If the PAMELA/ATIC data is indicative of dark matter annihila- 
tion, however, then the galactic annihilation rate must be boosted by 0(1000). This large boost 
factor considerably enhances the potential for galactic gamma ray and neutrino signals from the 
dark sector. In this context, hints from the WMAP Haze may also point towards an annihilation 
explanation of PAMELA/ATIC. 

In this paper, we have explored the robustness of dark matter annihilation predictions by 
considering cascade scenarios where dark matter annihilates into new resonances that in turn 
decay in one or more steps into standard model leptons. These cascade annihilation scenar- 
ios are directly motivated by the PAMELA/ATIC data, since light resonances can enhance the 
galactic annihilation rate through nonperturbative effects and explain the lepton-richness of the 
annihilation through kinematic thresholds. 

We have shown that electron and muon cascades give reasonable fits to the PAMELA/ATIC 
data. As a rule of thumb, the best fit dark matter mass and boost factor both scale as 2 n for 
rz-step cascade decays. We then compared these best fit values to constraints from gamma rays 
and neutrinos. The gamma ray bounds from FSR can be weakened by an order of magnitude 
through cascade decays, although increasing the length of cascades does not further weaken the 
bounds. Neutrino bounds for dark matter annihilating into muons are robust to changing the 
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length of the cascade, which is particularly relevant for models with large branching fractions to 
muons such as the minimal axion portal. 

Assuming standard NFW or Einasto halo profiles, there is tension between a dark matter 
annihilation interpretation of PAMELA/ATIC and the non-observation of galactic gamma rays 
or neutrinos. Such tension does not invalidate a dark matter annihilation hypothesis since there 
is considerable uncertainty in the dark matter halo distribution and velocity profile, and the 
constraints are uniformly weaker for shallower halo profiles. For gamma rays in particular, the 
galactic center and galactic ridge constraints assume an understanding of the dark matter halo 
profile in the inner 100 pc of the galaxy, where there is considerable uncertainty. The dark matter 
halo profile in the inner 4" of the galaxy is even more uncertain, so we do not consider radio 
measurements of synchrotron to be constraining. Also, for both gamma rays and neutrinos, the 
bounds can be weakened if the astrophysical boost factor for electrons/positrons is larger than 
those for gamma rays and neutrinos. 

If a dark matter annihilation scenario is realized in nature with the boost factor suggested by 
PAMELA/ATIC, then one would expect future experiments to see a gamma ray or neutrino flux 
given standard halo assumptions. ANTARES |56j, IceCube [57J, and KM3NeT [SS] will greatly 
increase current sensitivity to upward-going muons resulting from galactic neutrinos. Future 
atmospheric Cerenkov telescopes as envisioned in Ref. [5H] will also improve the prospects of 
finding gamma rays from dark matter annihilation. While we did not include the effect of ICS in 
our gamma ray analysis, ICS is expected to be a dominant dark matter annihilation signal in the 
energy range available to the Fermi Gamma-ray Space Telescope |60j. Ultimately, one hopes that 
future experiments could probe the detailed energy spectra of dark matter annihilation products 
to distinguish between direct annihilation and the cascade scenarios considered here. 

Note added: While completing this paper, we became aware of Ref. [61] which considers similar 
issues. 
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A Cascade Energy Spectra 



In this appendix, we present formulae for the energy spectra used in the text. In general, the 
energy spectra of final state particles in cascade annihilations are functions of all the intermediate 
masses and helicities. In the limit of large mass hierarchies and scalar decays, however, the energy 
spectra greatly simplify, and we use these simplified formulae in our analysis. 

Consider cascading fields fa of mass (mi+i > 2m,j) and a final state if) with mass m^. 
Cascade annihilation occurs through fa + i — > fafa (i — 1,2, • • ■), and in the last stage, 4>\ decays 
into if) + X. Let the energy of if) in the 0i rest frame be E . Defining 

x = , e = (24) 

mi mi 

the if) energy spectrum is a function of xq and eo 

where eo < Xq < 1. In the case where dark matter x annihilates directly into if) + X, we can 
regard fa as the initial state of dark matter annihilation, XX- m this case dN^/dxo is the primary 
injection spectrum with mi = 2moM- 

Now consider the previous step in the cascade annihilation, fa —> <fiifa, with one of the (pi 
decaying into if> + X. Let the energy of ip in the fa rest frame be E\ and define 

2£i 2m, 

x x = , ei = . (26) 

m 2 m 2 

Assuming isotropic scalar decays, the if) energy spectrum in the fa rest frame is 

^ = f 1 d cos 6 tdx ^ 5 ( 2x x -x - cos 6yjxl - e 2 ^l - ef) , (27) 

where 9 is the angle between the if) momentum and the 0i boost axis as measured in the 0i rest 
frame. 

Equation (1271) is complicated to solve in general, but in the limit e; — > (i = 0, l,---), it 
reduces to a simple convolution: 

dN^ fidxo dN^ 2 



dxi J Xl x dx ^ * ' ^ ^ 

where < x\ < 1 up to 0[ef) effects. This convolution can be iterated as many times as necessary 
to build up the desired energy spectrum for an n-step cascade decay: 

dNj, _ f 1 dx n ^ dN^ + g ^ 



dx n J x n %n—l d^n—l 
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where x n _i = 2E n _i/m n with E n _i being the energy of if) in the <p n rest frame, and < x n < 1 
up to 0(ef) effects. Note that we here adopt the normalization convention of 

1 . dN^ 



dx n -^ = l, (30) 
o ax n 

regardless of the value of n, so that the injection spectra per dark matter annihilation must be 
multiplied by the multiplicity of ip in the final state. 

A.l Direct electron spectra 

Here we derive the spectra of electrons/positrons arising directly from (f>i decay, 4>i ~^ e + e~ (or 
dark matter annihilation, XX ~^ e + e~). Ignoring the effect of final state radiation to smooth the 
spectrum, the electron energy spectrum is given by 

' 1X 6{l-x Q ), (31) 



dx 



o 



where we have adopted the convention ^dx a 5{l — x ) = 1. The positron energy spectrum is 
identical. 

Applying the simplified convolution formula in Eq. ( 1291) for an n-step cascade annihilation, we 
then find 



dN e 
dx\ 



1, (32) 



dN e 1 

— - = In—, 33 

dx 2 x 2 

,/A ' = Qn(Xn), (34) 



1 / l\ n_1 

Q n (x) = - - In- . (35) 



where we have defined 

()( r\ = 

(n — 1)! V x 

Note that these are energy spectra for one of the electrons (or positrons), so that the electron (or 
positron) injection spectra per dark matter annihilation XX ~^ 20„ 2 rt (e + e~) are 

^ = 2"^, (36) 

where x n = £" e / m DM with E e being the electron (positron) energy in the center-of-mass frame for 
dark matter annihilation. The direct annihilation case, XX ~ y e + e~, corresponds to n — 0. 
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A. 2 Electron and neutrino spectra from muon decay 

Here we discuss the spectra of electrons, positrons and neutrinos arising from muon decay. Con- 
sider 0i — > fi + fi~ (or XX ~ ¥ ^ + T-O followed by fi — > ez/ e z/ M . (One of the neutrinos here should 
be an anti-neutrino. We omit the particle-antiparticle identification here and below.) Assuming 
a massless electron, the (unpolarized) spectra of electrons and neutrinos in the rest frame of the 
muon are 

= dJ T^ = 6(^ 1 ) 2 -4(x_ 1 ) 3 , (37) 
ax -i ax -i 

= 12(x_0 2 - I2(x^f, (38) 

ax -i 

where we are using the notation x_i = 2_E_ 1 /m At , with E_i being the energy in the muon rest 
frame. 

Applying the cascade convolution for electrons and muon neutrinos 

~ l-^ri m 

dN^e dN.^ Uu 65 3 2 4 „ 19 , , 5^ . N 

= ~dtr = To8 _ 4^ + 27^ _ rs Q ^ + ~ 3 Q ^ (41 > 

and for electron neutrinos 

AT 

(42) 
(43) 
(44) 

Again, these are energy spectra for one of the electrons, positrons or (anti-)neutrinos. To obtain 
the injection spectra per dark matter annihilation, we must multiply the multiplicity factor, 2 n 
for n-step, and set x n = E/ijibm- Here, E is the energy of a particle in the center-of-mass frame 
for dark matter annihilation. 

For comparison, the corresponding formulae in the approximation of isotropic three-body 
decays are 

dN 

a ^ e ~ 2x_i, (45) 



dN^ Ve 
dx 


= 2- 


Gxq ~\~ 4xq, 




dN^ Ve 
dxi 


5 
~3 


+ 3x 1 — —x 1 


+ 2Q 2 (x 1 ), 


dN^ Ue 

dX2 


19 
18 


3 2 4 3 
2 ^2 + g ^2 


- ^Q2(x 2 ) + 2Q 3 {x 2 ) 



dx_i 
dN, 



n+l 



I) n+I [2x n + 2J2(-1)UM) = Q n+1 (x n ). (46) 



dx n \ | 
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As a rough rule of thumb, the electron spectrum for an n-step muon cascade has a shape between 
(n + 1)- and (n + 2)-step electron cascades. 



A. 3 Gamma ray spectra from final state radiation 

Primary gamma rays come from final state radiation in the decay <p% —>■ £ + £~7. In principle, one 
could do an exact calculation to C(«em) of the gamma ray spectrum, which would have the full 
eo = 2mg/mi dependence. Since we are using the simplified convolution formula in Eq. (1291) . it 
is not consistent to keep 0{e^) corrections in the exact gamma ray calculation, and it suffices to 
use twice the Altarelli-Parisi splitting formula 

dN^ = « em 1 + (1-xq) 2 r + ^ / 4(l-x ) \ 1 

dx 71 X \ \ C / J ' 

where the normalization of iV 7 is such that $dx n dN 1 / dx n gives the (average) number of photons 
per 0i decay. Note that the expression of Eq. (T4T1) becomes negative at xo > 1 — eel/4, which 
does not correspond to the kinematic threshold. The error from this, however, is formally an 
0{el) effect. 

Applying the simplified convolution formula, we obtain 

-l + lniW(xi)+£i(a;i)l, (48) 



dN^ 


«EM 


1 


dxi 


71 


• r l 


dN^ 


«EM 


1 


dx 2 


71 


•'■2 



where 



-l + \n^)R 2 (x 2 ) + S 2 (x 2 )), (49) 



R^x) = 2- x - x 2 + 2x\nx, (50) 

R2{x) = 2 — 3x + x 2 + x In x — a; (In x) 2 , (51) 

Si (a;) = - l\x + x 2 + 2xlnx + (2 - x - x 2 ) ln(l - x) - 2xLi 2 (x), (52) 

vr 2 _ _ ^,^\ _ 9 /7r 2 , 



S 2 ( x ) = i^— + 2 - 2C(3)Jx - 2x 2 - (^y -3jxlnx+ (2-3x + x 2 )ln(l -x) 

— x(lnx) 2 — xLi 2 (x) + 2xLi 3 (x). (53) 

The photon injection spectra per dark matter annihilation are then given by 

<±N 2 = r dN x< (54) 

where x n = E^/m^M, with being the photon energy in the center-of-mass frame for dark 
matter annihilation. Here, eo = rni/m-DM for direct annihilation and eo = 2vrnjm\ otherwise. 
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For the hardest gamma rays near x n — > 1, the behavior of dN^/dXn is 

fl * ^llnf 4(1 - Xo) V (55) 

dX 7T X \ 6 / 

^7 _ <*em!-Si ln / 4(l-xQ \ ^ 
dx\ Ti X\ y Cq y 

diV 7 aEM(l-^2) 2 , /4(l-x 2 )\ , . 

"1 o 111 2 • 57 

dx 2 7T 2x 2 V e o / 

Compared to direct annihilation into leptons, a 1-step cascade annihilation gives a gamma ray 
spectrum that is suppressed not only by ln(2m^/m 1 )/ln(m£/mDM) but also by an additional 
suppression factor of (1 — x) for the highest energy gamma rays. 

A. 4 Gamma ray subtlety for muons 

There are actually two contributions to the gamma ray spectrum for <px — ► /i + /i~7. In addition to 
final state radiation from muons, there is the radiative decay of the muon \x — > ev e v^. Formally, 
this contribution is suppressed by a factor of 1/ ln(m M /mi) or (1 — x n ) 2 , but for mi ~ m M , it is 
an important effect. 

The gamma ray spectrum in the muon rest frame is known in the limit that r = m 2 /m 2 is 
small [52] ■ Assuming unpolarized muons, we can derive the 0-step cascade annihilation spectrum 
from the muon rest frame spectrum 

^ ~ ^ 1 ^_x(x_0 In I + tMs-O V (58) 



dx-\ 3n x_i V r 



(ixo 37r xq V r 



To(xo)ln- + [/ (x ) , (59) 



where 



T_i(x) = (1 -x)(3-2x + 4x 2 -2x 3 ), (60) 
2 2 

T (x) = 3 + —x — 6x 2 + 3x 3 — -x 4 + 5xlnx, (61) 

/ 17 93 101 ^5 \ 

U^x) = (1-x) -— + — x-— x 2 + — x 3 + (3-2x + 4x 2 -2x 3 )ln(l-x)), (62) 
V 2 6 12 12 / 

j j i \ 17 3 191 2 23 3 7 4 / 2 D 2 o 3 2 A 1 M . 

Unix) = xH x x H — x + 3H — x — 6x + 3x x mil — x) 

ov ' 2 2 12 3 4 V 3 3 / V ; 

28 

— — xlnx + 5xln(l — x) lnx + 5xLi 2 (l — x). (63) 
The convolutions for 1- and 2-step decays are straightforward to derive. 
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Note again that Jdx n dN^ 7 /dx n give the (average) number of photons per muon decay. 
The photon spectra from radiative muon decay per dark matter annihilation are then given by 
2 n dN fJ- ^ 1 /dx n . The total photon injection spectra per dark matter annihilation are given by 

dN y dN Ur -r 1 \ , \ 

irl + 2 ^tf)< < 64 > 

where x n = E^/m^M, with E 1 being the photon energy in the center-of-mass frame for dark 
matter annihilation. 
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A. 5 Rare modes in the axion portal 

In Section [61 we consider bounds on rare a — > 77 and a — > 7r + 7r~7r decay modes in axion 
portal models. The axion portal has both a 1-step and a 2-step component, and the gamma ray 
spectrum for each can be calculated straightforwardly. The a — > 77 spectra are identical (up to 
normalization) to the 4>i — * e + e~ spectra already calculated: 



(65) 



dN a ^ 
dx\ 

= 21b—, (66) 
dx 2 x 2 

where the normalization of iV a ^ 7 is such that fdx n dN a ^/dx n gives the number of photons per 
a decay. 

For a — ► 7T + 7r _ 7r° followed by 7r° — ► 77, we can use the Q n function from Eq. ( 1461) if we assume 
that m a ^> m n and that the a — > 3n decay is isotropic: 

dNa ^ * 2g 2 (x x ), (67) 



1 

dN y°^ ~ 2Q 3 (x 2 ), (68) 
axo 



where again the normalization of A^ r a ^ 7r o^ 7 is such that Jdx n dN a ^o_^ 7 /dx n gives the number of 
photons per a decay. One should keep in mind that m a ~ 3m, in the region of interest, but the 
hierarchical cascade approximation is still reasonably representative of the true energy spectrum. 



B Leptonic Axion Portal 

In the minimal axion portal construction, the axion has large hadronic couplings, and is therefore 
strongly constrained by beam dump and rare meson decay experiments. In particular, the axion 
is forced to decay primarily into muons, and, as we saw in Section [51 there is some degree of 
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tension between a muon annihilation scenario and the absence of galactic neutrinos. Also, we saw 
in Section [6] that there are strong gamma ray bounds on the a — > 7r + 7r~7r° decay mode, which in 
the minimal axion portal arises from axion-pion mixing. 

Since the decay properties of the axion are irrelevant for dark matter freezeout, we can easily 
modify the couplings of the axion to standard model fields without losing the good features of 
this scenario. In particular, we can construct a leptonic axion portal model where the axion has 
no hadronic couplings. While such a leptonic axion could decay into muons as in the minimal 
axion portal, in the text we consider the less constrained case where the leptonic axion decays 
primarily into electrons. 

The simplest example for the leptonic axion portal can be constructed as follows. Vector- 
like fermion dark matter tfj/ip c obtains a mass from spontaneous symmetry breaking through the 
vacuum expectation value of a complex scalar Si: 

C = -iS^r + h.c, St = (f- a + e ia ^% (69) 

where ag is the pseudoscalar axion, sg is a light scalar, and f a is the axion decay constant, which 
is assumed to be of order TeV. As in the minimal axion portal, the masses of sg and ag can be 
considered as free parameters. In order for ag to decay into leptons, Si must be charged under 
a leptonic symmetry (which is softly broken to give a mass to ag). This requires introducing 
separate electron-type and neutrino-type Higgses: 

C = -\ e ih e e c - \ v ih v v c - AiS e h e h v + h.c. (70) 

These interactions force the standard model leptons to carry axial leptonic charges. (Small neu- 
trino masses can be obtained through the standard see-saw mechanism, and it is straightforward 
to extend the model to incorporate supersymmetry.) 

In order to eliminate the hadronic couplings of ag, the standard model quarks are assumed 
not to carry charges under the leptonic symmetry. This requires introducing one or more new 
Higgses for the quark sector, which must also be singlets under the leptonic symmetry. As long 
as the potential for these Higgses preserves the leptonic symmetry, then ag will have no hadronic 
couplings and cannot mix with the neutral mesons after the leptonic symmetry is spontaneously 
broken. 

The absence of hadronic couplings allows m ae to be lighter than 2m M , and thus ag to decay 
dominantly into electrons while satisfying beam dump and rare meson decay constraints. Also, 
the leptonic axion ag has no ag — > 7r + 7r~7r° decay mode, reducing the gamma ray constraints in the 
case that m ae > 3m n (where ag dominantly decays into muons). Since the strongest astrophysical 
bounds on light degrees of freedom come from hadronic couplings [63J , m ae might even be as light 
as 2m e , although a detailed study of the constraints on the leptonic axion is beyond the scope of 
this paper. 
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As an example of the quark sector interactions, there could be separate up-type and down- 
type Higgses. In this case, it is natural to assume a hadronic symmetry under which the quarks 
carry axial charges: 

C = -\ u qh u u c - \ d qh d d c - A q S q h u h d + h.c. (71) 

The axion contained in the field S q could then be the QCD axion and solve the strong CP 
problem. To avoid astrophysical constraints, however, the vacuum expectation value of S q must 
be ~?Z 10 10 GeV, much larger than (Se) = O(TeV). Explicit breaking of the hadronic symmetry 
must also be much smaller than that of the leptonic symmetry. 
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